/*
 *
 * Copyright (c) Kresimir Fresl 2003
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * Author acknowledges the support of the Faculty of Civil Engineering,
 * University of Zagreb, Croatia.
 *
 */

/* for UMFPACK Copyright, License and Availability see umfpack_inc.hpp */


#ifndef BOOST_NUMERIC_BINDINGS_UMFPACK_OVERLOADS_HPP
#define BOOST_NUMERIC_BINDINGS_UMFPACK_OVERLOADS_HPP

#include <Core/Utils/numeric/bindings/umfpack/umfpack_inc.hpp>
#include <Core/Utils/numeric/bindings/detail/array.hpp>
#include <Core/Utils/numeric/bindings/traits/detail/utils.hpp>
#include <cassert>

namespace boost
{
    namespace numeric
    {
        namespace bindings
        {
            namespace umfpack
            {
                namespace detail
                {
                    ////////////////////////////
                    // UMFPACK primary routines:
                    ////////////////////////////

                    // symbolic

                    inline
                    int symbolic(int n_row, int n_col,
                                 int const* Ap, int const* Ai, double const* Ax,
                                 void** Symbolic, double const* Control, double* Info)
                    {
                        return umfpack_di_symbolic(n_row, n_col, Ap, Ai, Ax,
                                                   Symbolic, Control, Info);
                    }

                    inline
                    int symbolic(int n_row, int n_col,
                                 int const* Ap, int const* Ai, traits::complex_d const* Ax,
                                 void** Symbolic, double const* Control, double* Info)
                    {
                        int nnz = Ap[n_col];
                        bindings::detail::array<double> Axr(nnz);
                        if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Axi(nnz);
                        if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Ax, Ax + nnz,
                                                    Axr.storage(), Axi.storage());
                        return umfpack_zi_symbolic(n_row, n_col, Ap, Ai,
                                                   Axr.storage(), Axi.storage(),
                                                   Symbolic, Control, Info);
                    }


                    // numeric

                    inline
                    int numeric(int, int,
                                int const* Ap, int const* Ai, double const* Ax,
                                void* Symbolic, void** Numeric,
                                double const* Control, double* Info)
                    {
                        return umfpack_di_numeric(Ap, Ai, Ax, Symbolic, Numeric, Control, Info);
                    }

                    inline
                    int numeric(int, int n_col,
                                int const* Ap, int const* Ai, traits::complex_d const* Ax,
                                void* Symbolic, void** Numeric,
                                double const* Control, double* Info)
                    {
                        int nnz = Ap[n_col];
                        bindings::detail::array<double> Axr(nnz);
                        if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Axi(nnz);
                        if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Ax, Ax + nnz,
                                                    Axr.storage(), Axi.storage());
                        return umfpack_zi_numeric(Ap, Ai, Axr.storage(), Axi.storage(),
                                                  Symbolic, Numeric, Control, Info);
                    }


                    // solve

                    inline
                    int solve(int sys, int,
                              int const* Ap, int const* Ai, double const* Ax,
                              double* X, double const* B, void* Numeric,
                              double const* Control, double* Info)
                    {
                        return umfpack_di_solve(sys, Ap, Ai, Ax, X, B, Numeric, Control, Info);
                    }

                    inline
                    int solve(int sys, int n,
                              int const* Ap, int const* Ai, traits::complex_d const* Ax,
                              traits::complex_d* X, traits::complex_d const* B,
                              void* Numeric, double const* Control, double* Info)
                    {
                        int nnz = Ap[n];
                        bindings::detail::array<double> Axr(nnz);
                        if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Axi(nnz);
                        if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Ax, Ax + nnz,
                                                    Axr.storage(), Axi.storage());
                        bindings::detail::array<double> Br(n);
                        if (!Br.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Bi(n);
                        if (!Bi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(B, B + n,
                                                    Br.storage(), Bi.storage());
                        bindings::detail::array<double> Xr(n);
                        if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Xi(n);
                        if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;

                        int status = umfpack_zi_solve(sys, Ap, Ai,
                                                      Axr.storage(), Axi.storage(),
                                                      Xr.storage(), Xi.storage(),
                                                      Br.storage(), Bi.storage(),
                                                      Numeric, Control, Info);
                        if (status != UMFPACK_OK) return status;
                        traits::detail::interlace(Xr.storage(), Xr.storage() + n,
                                                  Xi.storage(), X);
                        return status;
                    }


                    // free_symbolic

                    inline
                    void free_symbolic(double, int, void** Symbolic)
                    {
                        umfpack_di_free_symbolic(Symbolic);
                    }

                    inline
                    void free_symbolic(traits::complex_d const&, int, void** Symbolic)
                    {
                        umfpack_zi_free_symbolic(Symbolic);
                    }


                    // free_numeric

                    inline
                    void free_numeric(double, int, void** Numeric)
                    {
                        umfpack_di_free_numeric(Numeric);
                    }

                    inline
                    void free_numeric(traits::complex_d const&, int, void** Numeric)
                    {
                        umfpack_zi_free_numeric(Numeric);
                    }


                    ////////////////////////////////
                    // UMFPACK alternative routines:
                    ////////////////////////////////

                    // default control

                    inline
                    void defaults(double, int, double* Control)
                    {
                        umfpack_di_defaults(Control);
                    }

                    inline
                    void defaults(traits::complex_d const&, int, double* Control)
                    {
                        umfpack_zi_defaults(Control);
                    }


                    // symbolic with specified column preordering

                    inline
                    int qsymbolic(int n_row, int n_col,
                                  int const* Ap, int const* Ai, double const* Ax,
                                  int const* Qinit,
                                  void** Symbolic, double const* Control, double* Info)
                    {
                        return umfpack_di_qsymbolic(n_row, n_col, Ap, Ai, Ax, Qinit,
                                                    Symbolic, Control, Info);
                    }

                    inline
                    int qsymbolic(int n_row, int n_col,
                                  int const* Ap, int const* Ai, traits::complex_d const* Ax,
                                  int const* Qinit,
                                  void** Symbolic, double const* Control, double* Info)
                    {
                        int nnz = Ap[n_col];
                        bindings::detail::array<double> Axr(nnz);
                        if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Axi(nnz);
                        if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Ax, Ax + nnz,
                                                    Axr.storage(), Axi.storage());
                        return umfpack_zi_qsymbolic(n_row, n_col, Ap, Ai,
                                                    Axr.storage(), Axi.storage(),
                                                    Qinit, Symbolic, Control, Info);
                    }


                    ///////////////////////////////////////
                    // UMFPACK matrix manipulation routines
                    ///////////////////////////////////////

                    // triplet (coordinate) to compressed column

                    inline
                    int triplet_to_col(int n_row, int n_col, int nz,
                                       int const* Ti, int const* Tj, double const* Tx,
                                       int* Ap, int* Ai, double* Ax, int* Map)
                    {
                        return umfpack_di_triplet_to_col(n_row, n_col, nz,
                                                         Ti, Tj, Tx,
                                                         Ap, Ai, Ax, Map);
                    }

                    inline
                    int triplet_to_col(int n_row, int n_col, int nz,
                                       int const* Ti, int const* Tj,
                                       traits::complex_d const* Tx,
                                       int* Ap, int* Ai, traits::complex_d* Ax,
                                       int* Map)
                    {
                        assert(Tx == 0 && Ax == 0 || Tx != 0 && Ax != 0);
                        double *Txr = 0, *Txi = 0;
                        if (Tx != 0)
                        {
                            bindings::detail::array<double> ATxr(nz);
                            if (!ATxr.valid()) return UMFPACK_ERROR_out_of_memory;
                            bindings::detail::array<double> ATxi(nz);
                            if (!ATxi.valid()) return UMFPACK_ERROR_out_of_memory;
                            Txr = ATxr.storage();
                            Txi = ATxi.storage();
                            traits::detail::disentangle(Tx, Tx + nz, Txr, Txi);
                        }
                        double *Axr = 0, *Axi = 0;
                        if (Ax != 0)
                        {
                            bindings::detail::array<double> AAxr(nz);
                            if (!AAxr.valid()) return UMFPACK_ERROR_out_of_memory;
                            bindings::detail::array<double> AAxi(nz);
                            if (!AAxi.valid()) return UMFPACK_ERROR_out_of_memory;
                            Axr = AAxr.storage();
                            Axi = AAxi.storage();
                        }
                        int status;
                        status = umfpack_zi_triplet_to_col(n_row, n_col, nz,
                                                           Ti, Tj, Txr, Txi,
                                                           Ap, Ai, Axr, Axi, Map);
                        if (Ax != 0)
                        {
                            if (status != UMFPACK_OK) return status;
                            traits::detail::interlace(Axr, Axr + nz, Axi, Ax);
                        }
                        return status;
                    }


                    // scale

                    inline
                    int scale(int, double* X, double const* B, void* Numeric)
                    {
                        return umfpack_di_scale(X, B, Numeric);
                    }

                    inline
                    int scale(int n, traits::complex_d* X,
                              traits::complex_d const* B, void* Numeric)
                    {
                        bindings::detail::array<double> Br(n);
                        if (!Br.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Bi(n);
                        if (!Bi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(B, B + n,
                                                    Br.storage(), Bi.storage());
                        bindings::detail::array<double> Xr(n);
                        if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Xi(n);
                        if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;

                        int status = umfpack_zi_scale(Xr.storage(), Xi.storage(),
                                                      Br.storage(), Bi.storage(),
                                                      Numeric);
                        if (status != UMFPACK_OK) return status;
                        traits::detail::interlace(Xr.storage(), Xr.storage() + n,
                                                  Xi.storage(), X);
                        return status;
                    }


                    //////////////////////////////
                    // UMFPACK reporting routines:
                    //////////////////////////////

                    // report status

                    inline
                    void report_status(double, int, double const* Control, int status)
                    {
                        umfpack_di_report_status(Control, status);
                    }

                    inline
                    void report_status(traits::complex_d const&, int,
                                       double const* Control, int status)
                    {
                        umfpack_zi_report_status(Control, status);
                    }


                    // report control

                    inline
                    void report_control(double, int, double const* Control)
                    {
                        umfpack_di_report_control(Control);
                    }

                    inline
                    void
                    report_control(traits::complex_d const&, int, double const* Control)
                    {
                        umfpack_zi_report_control(Control);
                    }


                    // report info

                    inline
                    void report_info(double, int, double const* Control, double const* Info)
                    {
                        umfpack_di_report_info(Control, Info);
                    }

                    inline
                    void report_info(traits::complex_d const&, int,
                                     double const* Control, double const* Info)
                    {
                        umfpack_zi_report_info(Control, Info);
                    }


                    // report matrix

                    inline
                    int report_matrix(int n_row, int n_col,
                                      int const* Ap, int const* Ai, double const* Ax,
                                      int col_form, double const* Control)
                    {
                        return umfpack_di_report_matrix(n_row, n_col, Ap, Ai, Ax,
                                                        col_form, Control);
                    }

                    inline
                    int report_matrix(int n_row, int n_col,
                                      int const* Ap, int const* Ai,
                                      traits::complex_d const* Ax,
                                      int col_form, double const* Control)
                    {
                        int nnz = Ap[n_col];
                        bindings::detail::array<double> Axr(nnz);
                        if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Axi(nnz);
                        if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Ax, Ax + nnz,
                                                    Axr.storage(), Axi.storage());
                        return umfpack_zi_report_matrix(n_row, n_col, Ap, Ai,
                                                        Axr.storage(), Axi.storage(),
                                                        col_form, Control);
                    }


                    // report triplet (coordinate)

                    inline
                    int report_triplet(int n_row, int n_col, int nz,
                                       int const* Ti, int const* Tj, double const* Tx,
                                       double const* Control)
                    {
                        return umfpack_di_report_triplet(n_row, n_col, nz, Ti, Tj, Tx, Control);
                    }

                    inline
                    int report_triplet(int n_row, int n_col, int nz,
                                       int const* Ti, int const* Tj,
                                       traits::complex_d const* Tx,
                                       double const* Control)
                    {
                        bindings::detail::array<double> Txr(nz);
                        if (!Txr.valid()) return UMFPACK_ERROR_out_of_memory;
                        bindings::detail::array<double> Txi(nz);
                        if (!Txi.valid()) return UMFPACK_ERROR_out_of_memory;
                        traits::detail::disentangle(Tx, Tx + nz,
                                                    Txr.storage(), Txi.storage());
                        return umfpack_zi_report_triplet(n_row, n_col, nz, Ti, Tj,
                                                         Txr.storage(), Txi.storage(),
                                                         Control);
                    }


                    // report vector

                    inline
                    int report_vector(int n, double const* X, double const* Control)
                    {
                        return umfpack_di_report_vector(n, X, Control);
                    }

                    inline
                    int report_vector(int n, traits::complex_d const* X,
                                      double const* Control)
                    {
#if 0
      // see UMFPACK v 4.1 User Guide
      bindings::detail::array<double> Xr (n);
      if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
      bindings::detail::array<double> Xi (n);
      if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;
      traits::detail::disentangle (X, X+n,
                                   Xr.storage(), Xi.storage());
      return umfpack_zi_report_vector (n, Xr.storage(), Xi.storage(), Control);
#endif
                        return umfpack_zi_report_vector(n,
                                                        reinterpret_cast<double const*>(X),
                                                        reinterpret_cast<double const*>(0),
                                                        Control);
                    }


                    // report numeric

                    inline
                    int report_numeric(double, int, void* Numeric, double const* Control)
                    {
                        return umfpack_di_report_numeric(Numeric, Control);
                    }

                    inline
                    int report_numeric(traits::complex_d const&, int,
                                       void* Numeric, double const* Control)
                    {
                        return umfpack_zi_report_numeric(Numeric, Control);
                    }


                    // report symbolic

                    inline
                    int report_symbolic(double, int, void* Symbolic, double const* Control)
                    {
                        return umfpack_di_report_symbolic(Symbolic, Control);
                    }

                    inline
                    int report_symbolic(traits::complex_d const&, int,
                                        void* Symbolic, double const* Control)
                    {
                        return umfpack_zi_report_symbolic(Symbolic, Control);
                    }


                    // report permutation vector

                    inline
                    int report_perm(double, int, int np,
                                    int const* Perm, double const* Control)
                    {
                        return umfpack_di_report_perm(np, Perm, Control);
                    }

                    inline
                    int report_perm(traits::complex_d const&, int, int np,
                                    int const* Perm, double const* Control)
                    {
                        return umfpack_zi_report_perm(np, Perm, Control);
                    }
                }
            }
        }
    }
}

#endif // BOOST_NUMERIC_BINDINGS_UMFPACK_OVERLOADS_HPP
